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I. INTRODUCTION 

Solid-state materials deviate in various ways from the 
periodic idealizations sometimes used to describe them 
theoretically. In crystals, for example, such deviations 
can occur due to the presence of lattice distortions, im¬ 
purity atoms (that may or may not be magnetic), and va¬ 
cancies in the lattice, collectively termed as disorder. Dis¬ 
order can significantly influence material properties. A 
dramatic example for noninteracting electrons is Ander¬ 
son localization^ In spin systems, the focus of this work, 
quenched disorder in the spin interactions leads to frus¬ 
tration and can generate spin glasses^ A spin glass is a 
remarkable state of matter where, loosely speaking, spins 
are “frozen” in an irregular pattern, i.e., they display a 
very slow dynamics under external driving. Although 
this phase does not exhibit spin order in a traditional 
sense (e.g., as a ferromagnet), it is still distinctly differ¬ 
ent from a paramagnet. Experimentally, spin glasses are 
generally associated with a cusp in the ac susceptibility 
at a certain temperature above which the system behaves 
like a paramagnet^ Whereas no standard order param¬ 
eter captures a glassy transition, glass order parameters 
exist that do so (see, e.g., Refs. 0 and 0 ). 

In two-dimensional (2D) lattices (the focus of this 
work), the effects of quenched disorder on classical 
spins have been extensively studied over the years m 
However, the calculation of thermodynamic properties 
and correlation functions continues to be a computa¬ 
tional challenge While it is generally believed that no 
spin-glass phase exists for nonzero temperatures^ zero- 
temperature phases continue to be debated for various 
models of interest (see, e.g., Refs. 0 and [l). For quantum 
spin models, the sign problem^ in quantum Monte Carlo 
simulations and the exponential growth of the Hilbert 
space, relevant to full exact diagonalization calculations, 
represent an even greater challenge. Because of this, the 
properties of disordered quantum spin systems, and of 
quantum spin glasses in particular, have remained essen¬ 
tially unexplored. The existing literature in the subject 
has almost exclusively dealt with classical models. 


Our goal in this work is to use a recently introduced nu¬ 
merical linked cluster expansion (NLCE) for disordered 
systems^ to study the thermodynamic properties of the 
classical Ising (with S = 1/2), and quantum (spin-1/2) 
XY and Heisenberg models with bimodal random-bond 
disorder on the square and honeycomb lattices. NL- 
CEs allow us to obtain finite temperature properties of 
those models in the thermodynamic limit through the 
exact diagonalization of finite-size clusters. We specifi¬ 
cally study the energy, entropy, specific heat, and uniform 
magnetic susceptibility (for the magnetization in the z- 
direction) as a function of temperature. Since any glassy 
phase is only expected to emerge at zero temperature in 
these models, if at all, we do not study the spin-glass 
order parameter. In future work, we will study quan¬ 
tum quenches^— to examine the possibility of disorder 
driven localization in 2D. 

Our results are briefly as follows: For clean systems, 
they match well with known results for the square lattice, 
while we report additional results for the honeycomb lat¬ 
tice. For the disordered systems, we unveil some inter¬ 
esting features. In the Ising model, the uniform suscepti¬ 
bility x ~ 1 /T for all orders of the linked cluster expan¬ 
sion - we demonstrate this explicitly. The susceptibility 
in the Heisenberg model also increases with decreasing 
temperature (up to the lowest temperature we can ac¬ 
cess). These two cases differ starkly from the XY model 
where the uniform susceptibility (for magnetization in 
the z-direction, i.e., the same quantity calculated in the 
other two models) shows a plateau at low temperatures. 
At high temperatures, the clean and disordered models 
behave identically as regards the energy, specific heat, 
and entropy up to third order in inverse temperature, 
and show identical susceptibilities up to second order - 
we show this explicitly in Sec. HVI via a high temperature 
expansion. 

The presentation is organized as follows. In section 
m we introduce the three models we study (the spin- 
1/2 Ising, XY, and Heisenberg) and summarize some of 
their known properties in the square and honeycomb lat¬ 
tices. Section IIIII briefly describes NLCEs for systems 
with disorder. Numerical results for the latter, and their 
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comparison with those for clean systems, are presented 
in Sec. IIVI We conclude with a brief summary in Sec. [V] 

II. MODELS 

We are interested in thermodynamic properties of var¬ 
ious spin-1/2 models on the square and honeycomb lat¬ 
tices (see Fig. [[]). In the absence of disorder, and for 
nearest neighbor interactions, those models do not ex¬ 
hibit frustration on either lattice, which are both bipar¬ 
tite. While the thermodynamic properties of the vari¬ 
ous models studied here are qualitatively similar on both 
lattice geometries, there are significant quantitative dif¬ 
ferences, e.g., in the critical temperatures for the onset 
of quasi-long-range orders 14 These differences have their 
origin in the different coordination number in both lat¬ 
tices, with the honeycomb lattice having the smallest 
one. Hence, not surprisingly, for the spin-1/2 antifer¬ 
romagnetic Heisenberg model, the staggered magnetiza¬ 
tion is significantly suppressed on the honeycomb lattice 
as compared with the square latticed Disorder, on the 
other hand, leads to frustration on both lattices, and to a 
qualitative change of the intermediate to low temperature 
properties with respect to the clean systems. Frustration 
can be easily identified by trying to assign spins to the 
various sites in Fig. [T]to minimize the energy. If one takes 
the Ising Hamiltonian discussed below, one finds that for 
the overwhelming majority of disorder realizations there 
is no single spin configuration that minimizes the energy 
on all bonds. 

We should stress that both for the clean and disordered 
cases, we expect quantum fluctuations to strongly modify 
the results for the spin-1/2 XY and Heisenberg models 
from their classical counterparts (see, e.g., Ref. [l6j for 
examples of the effects of quantum fluctuations in frus¬ 
trated spin-1/2 XY and Heisenberg models on the hon¬ 
eycomb lattice). Our focus in this work is on systems 
with bimodal random-bond disorder, where the nearest- 
neighbor coupling between the spins takes values ± J with 
equal probability. 

A. Ising model 

The Hamiltonian for the spin-1/2 Ising model can be 
written as 

Rising = S*S] (1) 

(ij) 

where S'? = ±1/2 is the spin at site-i, and the sum is 
over nearest neighbors. In the absence of disorder (= 
J for all i,j ), Eq. {TJ has served as the quintessential 
model for magnetism and was solved exactly on a 2D 
square lattice by Onsager « 7 In the presence of continuous 
disorder, the Hamiltonian in Eq. CD, commonly known 
as the Edwards-Anderson model^ has also become a 
widely studied model for spin-glasses. 


The Ising model has a discrete Z 2 symmetry, i.e., the 
transformation S? —> —S'? for all j leaves the Hamil¬ 
tonian invariant. This symmetry can be spontaneously 
broken at sufficiently low temperatures to create an or¬ 
dered phase. For S? = ±1/2, the critical temperature is 
given by T C /\J\ = l/21og(l ± y/2) ~ 0.57 for the square 
lattice, and T C /\J\ = 1/2 log[(V3 ± lj/Cv^ — 1)] « 0.38 
for the honeycomb lattice, and is the same for the ferro- 
and antiferromagnetic cases*^ The latter is because, for 
bipartite lattices, a unitary transformation relates both 
models. This can be easily seen by rewriting the Hamil¬ 
tonian in Eq. 0 as Hi s i„g = T,i J i S i,A S i,B > where A 
and B are the two sub-lattice indices. One can then go 
from Ji — J ,, i.e., between the ferro- and antiferromag¬ 
netic models, via the transformation: S? A —> —S~ A and 

02 C2 

J i,B ^ °i,B' 

The specific heat of the clean model diverges at the 
critical temperature for both the square and the hon¬ 
eycomb lattices. Equivalently, the derivative of the en¬ 
ergy diverges at the critical temperature but the en¬ 
ergy remains finite throughout. For the antiferromag¬ 
netic model, the susceptibility is known to be finite ev¬ 
erywhere, with an infinite slope at the critical temper¬ 
ature. The maximum value of the susceptibility occurs 
at T m = 1.537 T c and T m = 1.688 T c for the square and 
honeycomb models respectively ; 14 : 19 i.e., above the crit¬ 
ical temperature. Our results for the clean systems are 
consistent with these. However, we cannot study the crit¬ 
ical phase or the properties of the system very close to 
criticality (see Figs. [3] and [4]). 

The Ising model with bimodal disorder has been exten¬ 
sively studied in the past 4— —— It is reasonably well 
established that no glassy phase exists for T > 0. 22 : 23 
It has, however, been established that a glassy phase 
appears at zero temperature in this model* 7 At finite 
temperature, our results for this case are described in 
Sec. Eva] 




FIG. 1. (Color online) Schematic of lattice models, square 
(left) and honeycomb (right), with bond disorder considered 
in this work. The red and blue bonds represent Jij = ±J. 
Black dots at vertices represent the spins. It is easy to see 
that the bond disorder causes frustration by trying to arrange 
the spins to minimize energy. 
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B. XY Model 

The spin-1/2 XY model can be written as 

i?x Y = Y, J ii&S? + S?SV) (2) 

(ij) 

where S^’ v are spin operators at site i, proportional to 
the Pauli matrices. We consider only the isotropic case, 
where the model has a continuous 17(1) rotation symme¬ 
try in the plane. The presence of a continuous symmetry 
precludes, via the Mermin-Wagner theorem^ a finite- 
temperature phase-transition involving the breaking of 
this continuous symmetry from occurring. For d < 2 di¬ 
mensions, the fluctuations in any putative ordered phase 
appearing from breaking a continuous symmetry grow 
with system size for finite temperatures, destroying any 
order (see e.g., Ref. |29|). Hence, this model has an or¬ 
dered phase only at T = 0. 

However, in 2D there can still be a finite temperature 
Berezinskii-Kosterlitz-Thouless (BKT) transitio n 30 ! 31 
below which the system exhibits quasi-long-range spin 
order. The critical temperature for the BKT transi¬ 
tion in the spin-1/2 XY model in the square lattice is 
T C /\J\ « 0.34i 32 i 33 We should stress that both classical^ 
and quantum models with continuous symmetries in 2D 
can exhibit this kind of behavior— We should also men¬ 
tion that most studies in the literature report results for 
the ferromagnetic XY model (J < 1). However, like for 
the Ising model, in the square and honeycomb lattice 
a unitary transformation relates the ferro- and antifer¬ 
romagnetic models and the critical temperature is the 
same in both. Our calculations for the susceptibility of 
the clean case on the square lattice converge down to tem¬ 
peratures of T/|J| ss 0.4, which is compatible with the 
onset of quasi-long-range order for T C /\J\ < 0.34i 32 i 33 

C. Heisenberg model 

The spin-1/2 Heisenberg model, also known as the 
XXX model, can be written as 

tfHeis = J ij Si ■ s j (3) 

<*j) 

where S, = (. Sf,Sf, Sf), and S^’ v,z are spin operators at 
site i, proportional to the Pauli matrices. The Heisen¬ 
berg model has an SU(2) symmetry, the highest of the 
three models considered in this work. The ground state 
in the clean case (Jij = J) is an ordered ferromagnet 
or antiferromagnet depending on the sign of the coupling 
constant J. Like for the XY model, long-range order only 
occurs at zero temperature. However, in contrast to the 
XY model, the 2D Heisenberg model does not develop 
quasi-long-range order at finite temperature. This is due 
to the fact that the internal symmetry group, SU( 2) 
[0(3)] for the quantum (classical) Heisenberg model, is 


non-Abelian, as opposed to the XY model which has 
an Abelian symmetry group, 17(1) [0(2)] for the quan¬ 
tum (classical) cases. Vortices or point-defects, which 
are responsible for the BKT transition^ occur only in 
the latter casej 36 Furthermore, in 2D, O(N) models with 
N > 3 or SU(N) models with N > 2, i.e., with non- 
Abelian symmetry groups, can be shown via perturbation 
theory to be asymptotically free, which for spin mod¬ 
els translates to a renormalization group flow towards 
paramagnetism. 37-38 

The results for the clean system presented here for the 
square lattice are nearly identical to the spin-1/2 results 
presented in Ref. [39], which also compares with other 
known results. 


III. NUMERICAL LINKED CLUSTER 
EXPANSIONS 

Numerical linked cluster expansions (NLCEs) are a 
computational technique that can be used to calcu¬ 
late extensive properties (per lattice site) of translation- 
ally invariant lattice systems. NLCEs, which are based 
on linked cluster expansions,—were introduced in 
Ref. [43|, where it was shown that the results obtained for 
thermodynamic properties were exact in the thermody¬ 
namic limit for systems with a finite correlation length. 
Furthermore, results could be obtained at significantly 
lower temperatures as compared to high-temperature ex¬ 
pansions for models that develop long-range order at zero 
temperature. In several subsequent works, NLCEs have 
been shown to be a powerful computational technique for 
determining not only thermodynamic properties of a va¬ 
riety of lattice models*^ - — but also for studying thermal- 
ization (or the lack thereof) at long times after a quench 
in isolated quantum systems—c— For completeness, we 
provide a brief description of NLCEs. Details of how to 
implement them can be found in Ref. [50l 

In NLCEs, the expectation value of an extensive ob¬ 
servable, per-site, O in a translationally invariant system 
can be calculated as a sum over contributions from all 
clusters c of different sizes that can be embedded on the 
lattice 

C7 = 55m(c)x Wo(c), (4) 

C 

where M (c) is a combinatorial factor equal to the number 
of ways that a particular cluster c can be embedded on 
that lattice. Wo(c) is the weight of cluster c for the 
given observable, which is calculated via the inclusion- 
exclusion principle 

W 0 (c) = 0(c)~J2 W o(s). (5) 

sCc 

O(c) is the expectation value of the observable O on the 
specific cluster c. Within NLCEs, 0(c ) is calculated us¬ 
ing full exact diagonalization. The expansion is carried 
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FIG. 2. (Color online) Clusters up to the fourth order in the 
site based NLCE on the square lattice. The two 3-site clusters 
have the same Hamiltonian. At fourth order, in addition to 
three clusters with identical Hamiltonian, two topologically 
new clusters appear - the closed loop and the “_L”. Each 
topologically new cluster is diagonalized separately. 


out order by order, i.e., by first considering clusters with 
one-site, then two-sites, and so on. Beyond the bare sum 
in Eq. ©, several resummation schemes exist that ac¬ 
celerate the convergence of the expansion— Here we will 
report results from Wynn and Euler resummation tech¬ 
niques whenever they allow us to extend the convergence 
of the results to lower temperatures.- 44 

Examples of clusters up to fourth order in the site- 
based NLCE used here on the square lattice are shown 
in Fig. [2] At third order, although there are two geomet¬ 
rically different clusters, they are topologically identical. 
They have the same Hamiltonian for the models with 
nearest interactions we consider here. At fourth order, 
there are three clusters (including the one with the four 
sites on a line) that again have the same Hamiltonian 
for the models considered here. However, two topologi¬ 
cally new clusters appear, namely, the closed loop and 
the “_L”. They have to be individually diagonalized. 
From the fourth order and beyond the number of distinct 
topological clusters increases rapidly (exponentially with 
the number of sites) , making the calculations increasingly 
costly. References 0 and |49| provide details on the var¬ 
ious topological clusters on the square and honeycomb 
lattices, respectively, as well as the number of such clus¬ 
ters as a function of the order of the expansion. 

Recently, in Ref. El it was shown that NLCEs can also 
be used to study systems with disorder. As described 
above, NLCEs can only be used for translationally in¬ 
variant systems. A priori , disorder breaks translational 
invariance. However, we are only interested in disorder 
averaged physical quantities. If we take a disorder av¬ 
erage over all possible disorder configurations in models 
with bimodal disorder, we restore translational invari¬ 
ance and the equivalent of Eq. © reads 

0 = Yl M{c) x Wo{c) (6) 


where the overline denotes the disorder average. The 


disorder average of the weights is in turn given by 

W^c) = 0(cj-£wb(5). (7) 

sGc 

In other words, the disorder average can be carried out 
order by order for each observable. The calculations then 
proceed as for the translationally invariant system if one 
replaces 0(c) by O(c). 

The computations in the presence of disorder are much 
more challenging than for translationally invariant sys¬ 
tems because of the additional average over all possible 
disorder realizations. For example, the largest clusters 
we consider here for the quantum models in the square 
lattice have 13 sites. At this order, there are a total of 
~ 1.9 x 10 6 connected clusters, of which 5,450 are topo¬ 
logically distinct— Each of these has to be fully diago¬ 
nalized for the 2 e disorder configurations corresponding 
to the £ bonds in the cluster. This has to, of course, be 
carried out for all lower orders as well, each with a dif¬ 
ferent set of topologically distinct clusters and disorder 
configurations. For the clean systems, we report results 
for cluster with up to 15 sites for the square lattice. In 
that case one has to diagonalize 42,192 topologically dis¬ 
tinct clusters with 15 sites. 

NLCE calculations fail to converge when correlations 
in the thermodynamic limit extend beyond the largest 
clusters considered. Therefore, NLCEs cannot be used 
to calculate observables in phases with long-range order 
unless one tailors the expansion to account for those— 
Since disorder usually shortens correlations at low tem¬ 
peratures, NLCEs are particularly useful to study quan¬ 
tum disordered systems, despite the increase of compu¬ 
tational cost because of the disorder average. It will be¬ 
come apparent when we discuss our results for the various 
thermodynamic properties of interest in this work, that 
NLCEs converge to lower temperatures in disordered sys¬ 
tems when compared to clean systems. As mentioned be¬ 
fore, quantum Monte Carlo simulations in the presence 
of disorder are severely limited by the sign problem. 


IV. RESULTS 


In this section, we discuss the results of our NLCE 
based study of the three spin models described in Sec. [TT] 
on the square and honeycomb lattices. In what follows, 
we set J = 1 as the energy scale. For each model and 
lattice geometry, we report the energy (E), entropy (S), 
specific heat ( C v ), and uniform susceptibility for the mag¬ 
netization along the 2 -axis (y) as a function of tempera¬ 
ture. These quantities are defined as 


</t) m - (hy 

N ’ v NT 2 


log Z E ((g*) 2 ) ~ (S z ) 2 

N + T ’ X NT 


where the overline denotes a disorder average, the an¬ 
gle brackets denote the thermal expectation value in the 
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grand-canonical ensemble (at zero chemical potential), N 
is the number of sites, and T is the temperature. As men¬ 
tioned earlier, the disorder average is carried out over all 
disorder configurations at each order in the NLCE. In all 
cases, the disorder averaged results are compared with 
the ones in clean systems. 

For all observables, we report bare NLCE results for 
the highest two orders of our site-based NLCE expansion, 
which are determined by the number of sites l in the 
largest clusters studied. Namely, we report the results 
from Eq. © when the contributions of all clusters with 
up to l sites are added, where l takes the two largest val¬ 
ues in our calculations for each model in each lattice ge¬ 
ometry. We also report results using two different resum¬ 
mation schemes, indicated as Wynn„ and Euler n . The 
subscript n denotes the order of the resummation pro¬ 
cess (see Ref. for details). The resununation schemes 
allow us to access significantly lower temperatures than 
the bare results in some cases as indicated below. 

Before discussing each model and observable in detail, 
we review a few general observations for completeness 
and pedagogy. In all models and observables discussed 
here, the numerical results at intermediate to high tem¬ 
peratures in the presence of disorder are close to those of 
the corresponding clean system. Whereas this is obvious 
for temperatures so high that the first order correction 
to the infinite temperature result is negligibly small, we 
notice from our results in Figs. ©[ 8 ] that the observables 
in clean and disordered models are indistinguishable for 
temperatures as low as T = 2 (barring the susceptibility). 

In order to show why this is so, we expand the par¬ 
tition function for small inverse temperature (3 = 1/T, 
Z = Tr(e~ p& ) « Tr(l - (3H + (3 2 H 2 /2 +...). The mod¬ 
els we consider have only nearest-neighbor coupling, i.e., 
H = + Sj)/2], where we have in¬ 

cluded a magnetic field h as a source to calculate the 
uniform susceptibility in the z-direction. The most gen¬ 
eral two-site Hamiltonian that describes all models of in¬ 
terest here is given by H ^ = 7 (Sf/SJ + S^S^) + AS?Sj, 
which becomes the Ising model for 7 = 0,A = 1, the 
XY model for 7 = 1, A = 0, and the Heisenberg model 
for 7 = A = 1. With this in mind, to first order in 
/3, the high temperature expansion for Z can be writ¬ 
ten as Z = 2 N - /3£ , ; > IV ; + h(Sf + 5?)/2], 
where N is the number of lattice sites. Note first that 
Tr(5'f) = 0 (the Pauli matrices are traceless), and second, 
that Tr(Hij) = 0, so that the linear correction vanishes. 
Therefore, to first order in /3, the clean and the disor¬ 
dered system have the same partition function. This is 
true regardless of the type of disorder. 

To second order, after expanding H 2 , we have 


Z = 2 n + — Tr 


'y ( Jij JklHij Hj 
L (■ ij),(kl > 


kl~ 


+ h 2 Y J s z l s z ] +hY J JiAjSk 

i,j ( ij),k 


( 9 ) 


Let us treat the above terms one by one. In the first term, 
for (ij) ^ (kl), the trace is identically zero as shown 
above. For the case when say i ^ Z, but j = k, we 
effectively have a new Hamiltonian for three neighboring 
spins. It is easy to verify explicitly that this trace also 
vanishes. The only possibility left for a nonvanishing 
contribution is (ij) = (kl). The trace of second term 
in the brackets is nonzero only for i = j. In the third 
term, again for k ^ i or j. the trace is zero. For, say j = 
k , considering only the diagonal elements of the matrix, 
S?SjSj oc Sf, we see that the trace vanishes. Taking 
these into account, and writing log Z to 0(/3 2 ), we have 


log Z 
N 


= log 2 + 


/ 3 2 


2 • 2 N N 


Tr 


(ij) 


Nh 2 


( 10 ) 


For both, the clean system and the system with bimodal 
disorder, J? = J 2 . We therefore get, 


\°gZ B 2 

— =log2 + T 


J z TV 2 Hi 


2 


0(/3 3 ) (11) 


where Hi is the Hamiltonian for a two-site system and 
the subscript “ 2 ” on the trace indicates a trace over the 
Fock space of a two-site system. The disorder average 
does not change the above expression, and therefore to 
this order, the clean and the disordered systems behave 
identically. The energy, for instance, is given to this or¬ 
der by E = -(9 log Z/d(3)/N = —(3J 2 Tr 2 (IL|)/2, the 
specific heat is C v = j3 2 J 2 Tr 2 (iL|)/2, and the uniform 
susceptibility is y = /3/4. 

In fact, it is straightforward to check that the par¬ 
tition function in the clean and disordered systems re¬ 
mains the same at O(0 3 ), except for terms proportional 
to h 2 . In other words, the energy, entropy, and specific 
heat are the same for clean and disordered systems up 
to third order in /3, but the uniform susceptibility devi¬ 
ates. The following expression can be derived along the 
lines of Eqs. © m for the third order correction to the 
partition function: 


log Z 
N 


log 2 


J 2 Tr 2 H. 


2 -n 2 


/3 3 h 2 

12 


h 2 

f T 

JijTrzHi&S?. ( 12 ) 


There is no sum over i,j, which represent the two sites 
in a 2-site system. For the clean model one just needs to 
replace with J in the above formula, while for the dis¬ 
ordered one, the disorder average produces two terms for 
± J (which cancel each other, implying that disorder ex¬ 
tends the paramagnetic behavior in the susceptibility to 
lower temperatures). One can the see that for h = 0, all 
thermodynamic quantities studied here, except the sus¬ 
ceptibility, are identical in the clean and disordered cases 
up to 0(/3 3 ). One can further verify that this changes at 
fourth order, where differences emerge between clean and 
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disordered systems. An example of a term that makes a 
difference at fourth order is a square loop with four sites 
and four bonds (see Fig. 0. Even in the Ising case, the 
Hamiltonian H 12 H 23 H 34 H 41 = 1/64 has a nonzero trace 
with a product of four different Jij . 

We should stress that, in all models studied here in the 
presence of disorder, we find that there is a significant 
amount of residual entropy (when comparing with the 
clean systems) at the lowest temperatures we are able 
to access with NLCEs. This is a clear indication of the 
lack of order at those temperatures. The behavior of the 
entropy, coupled with a saturation of the energy observed 
at the lowest temperatures accessible to us, confirms that 
there are many energy levels close to each other at low 
energies. This is the hallmark of frustration. 


A. Ising model 

Figures [3] and 0] show the energy (a), entropy (b), spe¬ 
cific heat (c), and uniform susceptibility (d) for the dis¬ 
ordered spin-1/2 Ising model on the square and honey¬ 
comb lattices, respectively. For the square lattice, we 
also plot the exact results for E, S, and C v for the clean 
model. 17 ’ 52 Several approximate analytic estimates exist 
for the susceptibility, but there is no closed form expres¬ 



FIG. 3. (Color online) Spin-1/2 Ising model on the 
square lattice. Panels (a)-(d) show the energy, entropy, spe¬ 
cific heat, and uniform susceptibility vs T, respectively. Solid 
lines depict disorder averaged quantities, while dashed lines 
depict results for the clean system. Thin lines report bare 
results for the last two orders of the NLCE, while thick lines 
report the results of two resummation techniques. A thin con¬ 
tinuous line following the results of the resummations reports 
results for a lower order of the same resummation technique 
and is used to gauge their stability. The dotted vertical line 
marks the position of the phase transition. The dashed dot¬ 
ted line shows exact analytic results for the clean system in 
the thermodynamic limit. 


sion for all temperatures. 

Figures 0Ka) andBJa) show that, as mentioned before, 
a generic feature in the presence of disorder is that the 
energy tends to plateau at the lowest temperatures ac¬ 
cessible to us. In that regime, the entropy is significantly 
higher than in the clean systems [Figs. (3Kb) and 0Kb)]. 
Distinct to the Ising models, the sharp divergence in 
the specific heat in the clean case [see Figs. HKc) and 
0Kc)], which indicates the phase transition, is replaced by 
what appears to be smooth peak in the presence of disor¬ 
der. The maximum of that peak appears at temperatures 
lower than the critical temperature in the clean case. At 
higher temperatures, Eq. m gives results that agree for 
E, C v . and S down to T « 1. We note that NLCE re¬ 
sults for the disordered model are well converged down 
to T ss 0.2 to 0.3, while the NLCE results for the clean 
model converge to temperatures that are close to T c , and 
agree with the analytic results in the disordered phase. 

Results for the uniform susceptibility in Figs.HKd) and 
0 Jd) show that this quantity behaves very differently in 
the clean and disordered systems. In the disordered case 
it exhibits a 1/T behavior at all temperatures, both on 
the square and honeycomb lattices. An order by order 
linked cluster expansion reveals that the only nonvanish¬ 
ing contribution to the susceptibility in the disordered 
case comes from the single-site system, and is trivially 



FIG. 4. (Color online) Spin-1/2 Ising model on the hon¬ 
eycomb lattice. Panels (a)-(d) show the energy, entropy, 
specific heat, and uniform susceptibility vs T, respectively. 
Solid lines depict disorder averaged quantities, while dashed 
lines depict results for the clean system. Thin lines report bare 
results for the last two orders of the NLCE, while thick lines 
report the results of two resummation techniques. A thin con¬ 
tinuous line following the results of the resummations reports 
results for a lower order of the same resummation technique 
and is used to gauge their stability. Resummation results are 
not presented for the clean case as they do not extend the 
convergence to lower temperatures. The dotted vertical line 
marks the position of the phase transition. 
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proportional to 1/T. All higher order contributions van- disorder averaged log Z for the clusters with one, two, 

ish. We show this for the first few orders of the NLCE and three sites shown in Fig. [2] 

on the square lattice. Below are the expressions for the 

I 


log Z (1) 
log Z& 
log Z {3) 


log ( 
log ( 


0 h 0 h s 

e 2 + e 2 


1 

2 L 
1 


pj 


0(J+*h) 


0{J-Ah) N 


+ 0^-0 


0h 3/9 h 0(J + h) 

e 2 + e 2 + e 2 


a (j-m 


1 r / ffh -0 h -0(.J+3h ) 

- |k>g (^ 2 e 2 + 2 e 2 + e 2 


+ 0 

— 0(J — 3h) 


+ e" 




(HJ-h) ' 


+ 0^-0 


(13) 


The uniform susceptibility can be obtained from y = 
0~ 1 (d 2 log Z)/(dh) 2 evaluated at h = 0. We get for the 
above three orders, 





(14) 


Already, one can see that no new contributions appear 
at the higher orders. To confirm this, we calculate the 
weights of the three clusters (see Ref. 0 for details about 


multiplicities, etc.): 

pjy(i) _ y (i) _ 0 _ 

W.W = x (2) - 2W& =0 ( 15 ) 

= x (3) - 2 - 3W W = 0. 

Indeed, only the single-site cluster contributes. One can 
check this at higher orders, and for the honeycomb lattice 
as well. This is not the case for the XY and Heisenberg 
models discussed below. 

Earlier studies for the Ising model with bimodal disor- 



FIG. 5. (Color online) Spin-1/2 XY model on the square 

lattice. Panels (a)-(d) show the energy, entropy, specific 
heat, and uniform susceptibility vs T, respectively. Solid lines 
depict disorder averaged quantities, while dashed lines depict 
results for the clean system. Thin lines report bare results 
for the last two orders of the NLCE, while thick lines report 
the results of two resummation techniques. A thin continuous 
line following the results of the resummations reports results 
for a lower order of the same resummation technique and is 
used to gauge their stability. The dotted vertical line marks 
the position of the phase transition. 



FIG. 6. (Color online) Spin-1/2 XY model on the hon¬ 
eycomb lattice. Panels (a)-(d) show the energy, entropy, 
specific heat, and uniform susceptibility vs T, respectively. 
Solid lines depict disorder averaged quantities, while dashed 
lines depict results for the clean system. Thin lines report bare 
results for the last two orders of the NLCE, while thick lines 
report the results of two resummation techniques. A thin con¬ 
tinuous line following the results of the resummations reports 
results for a lower order of the same resummation technique 
and is used to gauge their stability. Note that the results 
converge to similar temperatures as in the square lattice. 































der on the square lattice found a low-temperature scaling 
of the specific heat (i.e., the exponent a in a power law 
fit of the low-temperature specific heat, C v ~ T~ a ) that 
is different from the model with continuous disorder^ 
At even lower temperatures, a crossover in the scaling 
behavior of C v has been reported;- Unfortunately, our 
results do not converge at low enough temperatures to 
observe such power laws. However, for T > 0.3, our re¬ 
sults are consistent with those in other studies^ (since 
we consider S z = ±1/2, our temperatures are lower by a 
factor of four from those studies, which took S z = ±1). 


temperature for low temperatures shows that increasing 
temperature does not increase the disorder in the spin 
correlations in the ^-direction. 

We should add that the classical XY model has been 
studied in the presence of Gaussian-random dilution^ 
and bimodal dilution^ of ferromagnetic bonds. In these 
works, the BKT transition was seen to slowly disappear 
as the dilution was increased. Here we only have consid¬ 
ered the fully disordered case, i.e., an equal distribution 
of ferro- and antiferromagnetic bonds, so we do not ex¬ 
pect that any remnants of the BKT phase are present in 
our calculations. 


B. XY model 

Figures E] and [6] show results for the spin-1/2 XY model 
on the square and honeycomb lattice, respectively. The 
results for all quantities are well converged down to about 
T ft 0.1 to 0.2. Figures [H(a),(b) and |6](a),(b) show that 
the behavior of energy and the entropy is qualitatively 
similar to the one observed in the Ising model. However, 
the results for the XY model in the presence of disorder 
converge at lower temperature than those for the Ising 
model. For the XY model, the specific heat in the pres¬ 
ence of disorder exhibits a peak that is well resolved by 
our NLCE [Figs. Etc) andEJc)]. We note that the energy, 
entropy, and specific heat follow the second order result 
in Eq. m for T > 2 in the square lattice and T > 0.7 
in the honeycomb lattice. 

Interestingly, Figs.Etd) andEtd) show that in the XY 
model the uniform (z-)susceptibility in the presence of 
disorder exhibits a plateau for low temperatures. This 
is qualitatively different from the behavior observed for 
the Ising model. The fact that the response to an ex¬ 
ternal magnetic field in the ^-direction is independent of 


C. Heisenberg model 

Figures [7] and [5] show results for the spin-1/2 Heisen¬ 
berg model on the square and honeycomb lattices, re¬ 
spectively. In Figs. [7](a) andE^a), one can see that the 
plateau in the energy at low temperatures is the clear¬ 
est of all models studied in this work. The onset of this 
plateau occurs at T ft 0.2 for both models. Remarkably, 
in the honeycomb geometry, the energy in the presence 
of disorder converges down to T ft 0.02. The entropy 
[Figs. [71[b) andEKb)] behaves similarly to the entropy in 
the XY model, but also converges to very low temper¬ 
atures (T ft 0.03 to 0.04) in the honeycomb geometry. 
Like for the XY model, the specific heat exhibits a clear 
peak in the presence of disorder. The temperature at 
which the maximum of that peak occurs is very close 
(but slightly larger) to the one in the clean model. 

In contrast to the XY model, the uniform suscepti¬ 
bility in the presence of disorder increases rapidly with 
decreasing temperature at the lowest temperatures ac¬ 
cessible to us. The susceptibility therefore behaves qual- 



T T 



FIG. 7. (Color online) Spin-1/2 Heisenberg model on 
the square lattice. Same as Fig. [5] but for the spin-1/2 
Heisenberg model. 


FIG. 8. (Color online) Spin-1/2 Heisenberg model on 
the honeycomb lattice. Same as Fig. [6] but for the spin- 
1/2 Heisenberg model. 
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itativcly similar to the Ising model. It is worth emphasiz¬ 
ing that our results for all observables in the honeycomb 
lattice appear to converge at temperatures significantly 
below T = 0.1, and quite close to T = 0.01 for the energy, 
entropy, and uniform susceptibility. 

Our calculations for the Heisenberg model reach lower 
temperatures and access regimes beyond what has been 
possible with quantum Monte Carlo simulations. For ex¬ 
ample, Ref. [55| studied the case of diluting an antifer¬ 
romagnetic model with ferromagnetic bonds of varying 
concentration. The results presented there did not reach 
the equal probability J = ±1 case discussed here because 
of the sign problem. For lower concentrations of ferro¬ 
magnetic bonds (below 30%), the calculations were still 
limited to temperatures above T « 0.3. 

V. SUMMARY 

We have used numerical linked cluster expansions to 
study thermodynamic properties of spin models with bi- 
modal (± J) bond disorder. The results reported are in 
the thermodynamic limit at temperatures for which they 
are well converged. 

We have unveiled various interesting effects of disorder 


in spin-1/2 Ising, XY, and Heisenberg models. For all 
models we find that disorder leads to a saturation of the 
energy at the lowest temperatures accessible to us, in a 
regime where the entropy is higher than in clean systems. 
This makes apparent that there are many low lying en¬ 
ergy states. For the disordered classical Ising model, on 
both the square and the honeycomb lattice, the diver¬ 
gence of the specific heat in the clean case is replaced by 
a peak, and the uniform susceptibility follows an inverse 
temperature law for all temperatures in the presence of 
disorder. This was explicitly verified order by order. In 
the Heisenberg model, we also find that the susceptibil¬ 
ity increases with decreasing temperature for all temper¬ 
atures accessible to us in the presence of disorder. In the 
XY model, on the other hand, we find that the suscep¬ 
tibility exhibits a plateau at low temperatures. On both 
the XY and Heisenberg models, our NLCE calculations 
were able to resolve a peak in the specific heat. 
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